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We explore the nature of the quantum phase transition between a magnetically ordered state with 
collinear spin pattern and a gapless Z2 spin liquid in the Heisenberg-Kitaev model. We construct 
a slave particle mean field theory for the Heisenberg-Kitaev model in terms of complex fermionic 
spinons. It is shown that this theory, formulated in the appropriate basis, is capable of describing the 
Kitaev spin liquid as well as the transition between the gapless Z2 spin liquid and the so-called stripy 
antiferromagnet. Within our mean field theory, we find a discontinuous transition from the Z2 spin 
liquid to the stripy antiferromagnet. We argue that subtle spinon confinement effects, associated 
with the instability of gapped U(l) spin liquid in two spatial dimensions, play an important role at 
this transition. The possibility of an exotic continuous transition is briefly addressed. 

PACS numbers: 



I. INTRODUCTION 

Theoretical understanding of various quantum spin liq- 
uid phases and their possible realization in candidate ma- 
terials consist one of the most important and interesting 
questions in modern condensed matter physics. 1 3 Much 
impetus in this field was derived from the discovery, by A. 
Kitaev, of an exactly solvable spin- 1/2 Hamiltonian on 
the Honeycomb lattice with a spin liquid ground state. 4 
This Hamiltonian, known as the Kitaev model (see be- 
low), has been studied intensively and it has been shown 
to support a gapped and a gapless Z2 spin liquid with 
fractionalized excitations. While the gapped phase has 
abelian anyon excitations, the gapless phase, in the pre- 
sense of appropriate perturbations, supports non-abelian 
anyons. 4 7 

Interestingly, it has been shown recently by Chaloupka 
et al. 8 that such a Kitaev model can indeed arise in 
layered Honeycomb lattice materials in the presence of 
strong spin-orbit coupling. In particular, they showed 
that in certain iridate magnetic insulators, the low en- 
ergy Hamiltonian for the pseudospin J = 1/2 iridium 
moments is given by a linear combination of the antifer- 
romagnetic Heisenberg model (Hn) and the Kitaev model 
(#k): 

H = (1 - a)H H - 2aH K (1) 

where a, expressed in terms of the microscopic parame- 
ters, determines the relative strength of the Heisenberg 
and the Kitaev interactions (the detailed forms of 
and Hk are given below). 

Subsequent to this suggestion, two honeycomb lattice 
compounds, Na2lr03 9 13 and L^IrOs 11 , have been dis- 
covered where such a Heisenberg-Kitaev (HK) model has 
been suggested to capture the low energy magnetic be- 
haviour. Meanwhile there have been intense numeri- 
cal studies 11,14,15 determining the phase diagram for the 
model as a function of a, magnetic field, including fur- 



ther neighbour interactions 16 and even doping 17,18 (we 
shall not consider the effect of the further neighbour in- 
teractions or doping in this paper). These studies reveal 
a rich phase diagram as a function of a which is shown 
in figure 1(a). There are three phases for a G [0,1] 8 . (1) 
The Neel Phase: At a = we have the antiferromagnetic 
Heisenberg model, which gives rise to collinear Neel order 
on the bipartite honeycomb lattice (figure 1(b)). As a is 
increased, the Neel state becomes unstable at a « 0.4 to 
a (2) Stripy order (figure 1(c)). The stripy state can be 
seen as antiferromagnetically coupled chains which are 
then coupled ferromagnetically. Between a « 0.4 — 0.8, 
the stripy phase is stable. Beyond a ~ 0.8 it gives way 
to a (3) gapless Z 2 spin liquid which is continuously con- 
nected to the gapless phase of the Kitaev model (for 
a = 1). While the phase transition between the Neel and 
the Stripy phase appears to be discontinuous, numerical 
studies including density matrix renormalization group 
(DMRG) 15 and exact diagonalization (ED) 8 results sug- 
gest that the transition between the spin liquid and the 
stripy state is continuous or weakly first-order. DMRG 
also indicates that turning on a magnetic field at the crit- 
ical point between the spin liquid and the stripy phase 
immediately opens up a polarized phase 15 and hence sug- 
gests that the phase transition between the spin liquid 
and the stripy phase may actually be governed by a multi- 
critical point. 

In this paper, we explore the nature of the phase tran- 
sition between the Z2 spin liquid and the stripy ordered 
phase as a function of a. Contrary to the original de- 
scription of the spins in terms of Majorana fermions 
employed by Kitaev, 4 we utilize a more conventional 
slave-particle approach to describe the Kitaev spin liq- 
uid which is then easily extended to include the Heisen- 
berg term. This helps us to describe the transition be- 
tween the stripy state and the spin liquid state within a 
slave-particle mean field theory. Our slave-particle for- 
mulation differs from that of Burnell et al. 19 and You et 
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FIG. 1: (a) The Phase diagram for the Heisenberg-Kitaev 
model and the spin pattern in the (b) Neel and (c) Stripy 
phases. 



by introducing the HK model in detail in Section II, and 
describe a change of basis 8 which we will use throughout 
the paper. This basis change allows us to capture the 
transition between the stripy magnet and the spin liquid 
easily. We then formulate the slave particle description 
of the model which we use to gain insight into the Kitaev 
model in Section III. We examine the exactly solvable Ki- 
taev limit of this model within this formulation, and show 
that the properties of the model are reproduced within 
our formulation. An examination of the gauge structure 
of the model follows, which we use to argue that the Z 2 
spin liquid phase is stable prior to the formation of mag- 
netic order. In the presence of magnetic order, however, 
the invariant gauge group is changed into a C/(l) struc- 
ture. In Section IV, we describe the mean field results 
in detail. We also argue that these results indicate that 
the transition is driven by confinement of the spinons, 
once we go beyond mean field theory. Finally, in Sec- 
tion V, we conclude with a discussion of our results and 
indicate the possibility of continuous transition induced 
by quantum fluctuations. The details of mean-field self- 
consistent equations are given in Appendix A. 



a/. 17 , allowing us to extend our analysis into the mag- 
netically ordered region by including a direct magnetic 
decoupling channel. Within our mean- field treatment, 
the transition appears to be first order with a discon- 
tinuous jump in the magnetic order parameter that is 
greater than that predicted by numerical calculations. 15 
However, we find that this transition is brought about 
by subtle non-perturbative effects associated with the 
confinement of a gapped U(l) spin liquid in two spatial 
dimensions. 20 In particular, we find, within mean-field 
theory, that on decreasing a from 1 the gapless Z 2 spin 
liquid goes into a gapped U(l) spin liquid phase with 
the simultaneous onset of magnetic order, albeit discon- 
tinuously. However, such a gapped U(l) spin liquid is 
unstable to non-perturbative instanton effects in two spa- 
tial dimensions, 20 which leads to immediate confinement 
of the spinons resulting in a conventional stripy order. 
We discuss the possible limitations of the present mean 
field theory and point out that non-perturbative quan- 
tum fluctuations beyond mean-field may allow a more 
exotic continuous transition. 

While, at present, it is not clear whether the HK model 
is relevant to the physics of the two above compounds 
in particular 21 (and in any case these compounds or- 
der magnetically), it presents an interesting microscopic 
Hamiltonian in several aspects. In addition to having a 
rich phase diagram of magnetically ordered and spin liq- 
uid phases, as discussed below, it offers an opportunity to 
study the effects of perturbation around an exactly solv- 
able Z 2 spin liquid (obtained at a = i) 15 > 22 > 23 . This latter 
direction allows us to study regular slave-particle mean- 
field theories 1 ' 16 19 ' 22, 24 in a more controlled setting. In 
fact, the slave-particle mean- field theory is expected to 
be exact at the exactly solvable point a = 1. 

The rest of the paper is organized as follows. We begin 



II. THE HEISENBERG-KITAEV 
HAMILTONIAN 

We start with the discussion of the Heisenberg-Kitaev 
Hamiltonian. The HK Hamiltonian is given by Eq 1, 
where the Heisenberg and the Kitaev 4 terms are given 

by 

ff H = ££-4 (2) 

(ij) 

H *= E E s " s r ( 3 ) 

(3=x,y,z {ij) ,(3 — links 

Si denotes spin 1/2 operators denned on the sites of the 
Honeycomb lattice, and (ij) denotes the nearest neigh- 
bour bonds. The Heisenberg term (Hn) is the usual spin 
rotation invariant antiferromagnetic Heisenberg Hamil- 
tonian, coupling spins on all nearest neighbour bonds. 
In contrast, the Kitaev term 4 (i^x) couples the x com- 
ponents of the spins on one of the directions of bonds 
(referred to as x — links) on the honeycomb lattice, the 
y components of spins on the y — links, and the z com- 
ponents on the z — links, as shown in Figure 2. More 
precisely, the Kitaev model that we have written down 
is the isotropic Kitaev model where the couplings on the 
x, y and z links are equal. 4 

The HK model does not have continuous spin rotation 
symmetry other than at the points a = and 0.5. This 
stems from the Kitaev part of the Hamiltonian which 
is devoid of continuous spin rotation symmetry. At this 
point it is useful to note the important symmetries of the 
HK model, which are 

1. ^ spin rotation about [111] spin axis along with 
Cs lattice rotations about any site. 
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FIG. 2: The different interactions of the Kitaev model 4 where 
the x,y and the z bonds are shown. The black and the 
white circles denote the two sublattices A and B. Rab — 
| [l, — and R V AB = \ [l, Vs] are the two unit vectors. 




FIG. 3: The rotated basis: the squares are left invariant, the 
circles are rotated about the z bonds, the triangles about the 
x bonds and the pentagons about the y bonds. This rotation 
was first described by Khaliulin 26 and Chaloupka et al 8 



2. Inversion about any plaquette center 

3. Inversion about any bond center. 

4. Time reversal. 

The Cs symmetry ensures that there are three differ- 
ent stripy phases, which we will refer to as the x, y and 
z stripy phases (see later). For the /3(= x,y,z) stripy 
phase, the spins are oriented along the j3 axis, with the 
(3 links being ordered ferromagnetically and the remain- 
ing two links ordered antiferromagnetically. Figure 1(c) 
shows one of the three possible stripy phases, namely z 
stripy phase. 

At the point a = 1, the model can be exactly solved 
by transforming the spins into products of Major ana 
fermions, with a background of frozen Z2 fluxes over 
plaquettes 4 . This is a gapless Z 2 spin liquid, with 
strictly nearest neighbour spin-spin correlations. 7 On the 
other hand, for a = we have the pure spin rota- 
tion invariant nearest neighbour Heisenberg antiferro- 
magnet where both numerical methods and semiclassical 
approaches give 2-sub-lattice Neel order. 25 In addition 
to these points, the model has another exactly solvable 
point at a — 0.5, where the stripy state is the exact 
ground state. 8 ' 26 This is easy to see by doing a selective 
rotation of the spins on the honeycomb lattice. It turns 
out that this rotated basis is useful to describe the transi- 
tion between the stripy phase and the spin liquid. Hence, 
we shall recall the the essence of the rotation as pointed 
out by Khaliulin 26 and Chaloupka et al 8 

A. The HK model in the rotated basis 

The transformation of the spin basis required to reveal 
the exactly solvable point at a = 0.5 is described in fig- 
ure 3. This transformation requires different spins to be 



rotated about different axis, depending on their position 
in the lattice as described in the figure. We first choose a 
set of spins which are positioned on third nearest neigh- 
bour sites at opposite corners of the hexagons throughout 
the lattice, and hold these spins fixed. We next rotate 
the three spins that are adjacent to these fixed spins by 
7r about the spin axis corresponding to the bond which 
connects it to the fixed spin. This has the net effect of 
transforming the Heisenberg term as 

H H -> -H H + 2H K (4) 

and leaving the Kitaev term invariant, i.e. 

H K ^H K . (5) 

From now on, we use spins in the rotated basis. However, 
for the sake of brevity, we shall continue to use the same 
symbol for the spins and the Hamiltonians. In this basis, 
the Hamiltonian (given by Eq. 1) becomes 

H->H=-(1- a)Hn - 4(a - ±)H K . (6) 

In this form, the exactly solvable point at a = 0.5 is 
clearly visible, as here the coefficient of the Hk term 
is zero and this is simply the ferromagnetic Heisenberg 
model with a ferromagnetic ground state in terms of 
the rotated spins. On undoing the rotations we recover 
the stripy-antiferromagnetic ordering in terms of the un- 
rotated spins. 8 ' 26 

Since we wish to particularly examine the transition 
between the Kitaev spin liquid and the stripy anti- 
ferromagnetic state, we find it easier to use this rotated 
basis. Also, it is helpful to think about deviations from 
the exactly solvable point at a = 0.5 in order to simplify 
the couplings in the region of interest. We achieve this 
by introducing the parameter 
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This gives 



H 



(8) 



Finally, we restrict ourselves to the region 5 G [0,0.5] 
where these couplings are purely ferromagnetic. We re- 
mind ourselves that S = now refers to the exact ferro- 
magnetic state (or the stripy state in the original basis) 
while 5 = 0.5 refers to the exactly solvable Kitaev point. 

We take this rotated model as the starting point of our 
slave-particle analysis. 



III. SLAVE PARTICLE FORMULATION 

Having written down the Hamiltonian (Eq. 8) in the 
desired form, we now introduce the slave-particle decom- 
position of the (rotated) spins. We write the spin- 1/2 
operator as a bilinear of two spin- 1/2 fermionic spinons 



as 



1,24 



s 1 ; = \fU^Uf^ 



(9) 



where fj a (& =t>i) are the fermionic spinon annihi- 
lation operators which satisfy regular fermionic anti- 
commutation relations. The above representation of the 
spin operators, along with the single fermion per site con- 
straint 



(10) 



constitutes a faithful representation of the spin- 1/2 
operators. 1,24 

The bilinear spin-spin interaction is a quartic term in 
the spinon operators. Within mean field theory, we now 
seek to decouple these quartic spinon terms into stable 
decoupling channels which are quadratic in terms of the 
spinons operators. In general, we need to keep both 
the particle-hole and particle-particle channels for the 
spinons. 1 

However, we note that both the terms in the final 
Hamiltonian (Eq. 8) have ferromagnetic interactions. 
Thus the usual decoupling 1 in terms of the spin-singlet 



particle-hole and particle-particle channels is unstable 
within a auxiliary field decoupling scheme. Instead, it 
was shown by Shindou and Momoi 27 that the correct 
spin liquid decoupling scheme for such interactions is into 
the spin triplet channels (both particle-hole and particle- 
particle). This is done as follows. We write the a-th 
component of the spin-spin interaction as 



(3=x,y,z 



e?Ie?+d?Id? 



l,p 1,1 



~ ~4 ' 
(") 

where S a ^ = 1(0) for a = /3(a ^ j3) (not to be confused 
with the parameter 5) is the Kronecker delta function, 



E i, P = ^fi+paWUpfip, 



D i,p = ~fi+pa[iv y ^]a(3fi(3 



(12) 
(13) 



and rii = 1 is the number of spinons per site. 

In addition to these hopping and pairing decouplings 
which capture the spin liquid, we introduce a direct chan- 
nel or magnetic decoupling, 



(14) 



which, without loss of generality, we choose to be in the 
S z direction. We include this decoupling explicitly in or- 
der to access the ferromagnetic state and due to the fact 
that it is the competing order in the spin liquid phase. It 
is important to note that when this operator has a non- 
zero expectation value (in the unrotated basis) it explic- 
itly breaks the the discrete symmetry corresponding to 
a lattice rotation by ^ about an individual site in con- 
junction with a spin rotation by 4^ about the [111] spin 
axis (refer to our discussion of the symmetries of the HK 
model). 

Using the above general ansatz, the mean-field spinon 
Hamiltonian for the rotated HK model is given by 
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MF 
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MF 
K •> 



(15) 



where 



K F = \ EE( m tfU ff *Wi./* + fl P J<T z Ufi +P ,e) - 2m 2 + {flJi, P ■ S a pf M + h.c.) - 2|£ iiP | 2 (16) 

% p 

+ (ft a Di, P ■ (-i™ v Uf! +p>0 + h.c.) - 2|A,p| 2 ), 

i p,r 
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Here i refers to the honeycomb lattice sites and p,r(= 
x,y,z) correspond to the link types and spin compo- 
nents respectively. Eq. 17 is the most general mean field 
Hamiltonian with the above decoupling channels. 



Taking the Fourier transform and restricting the pa- 
rameters E, D and m for each bond type to have the 
symmetry of the lattice, we get 



K F =2 ((fla t A e «f*Wfk t f},B + h.C.) + (fU A A a ,(k)f^ B + h.C.)) - ^ 5^(|4,| 2 + |4| 2 ) 



(18) 



m , 



k rj=A,B 



HM F =2 J2 ((fU A ^(k)f k ^ B + h.c.) + (/i a>A A aj9 (fc)/i fc|i9|jB Hh /i.e.)) - ^(1 - 8 p mk\ 2 + i^i 2 )> 



(19) 



r 



where iV s ^ e is the number of lattice sites and we have 
defined the Fourier transform of the spinons as: 



/i 



' Ka ' L ~ VN^ elk ' Rifi ' a ' L 



(20) 



(N is the number of unit cells, a =t, 4- and L = A, B is 
the sub-lattice index) and also introduced 



We can write the above mean field Hamiltonian (Eq. 
15) in a more compact way as a Bogoliubov-de-Gennes 
Hamiltonian for the spinons using the 4-component 
Nambu spinons, 



fi,t fi,l At fhi 



(27) 



= gma z a p, 

V 
p,r 



(21) 
(22) 

(23) 

(24) 



Aop(fe) = g I> " W^^I-urV]^, (25) 



where we denote the unit lattice vectors (refer to figure 
2) with 

Rab = (y ^B = (y' y)' R Z AB = ®- ( 26 ) 



The mean field Hamiltonian can now be written as 



1 $ 

(g - ^M/iU^W/i./S + /i+p,a[ f7 1«/9/i+P,/3)), 

(28) 



where 



C = %(3m 2 (i - <5) + - <5) + 25(1 - 5 p , r 



e;\ 2 + \d;\ 2 )) (29) 



and 



U ** = I E ( - ^ - *) - 4 ^ " V)) (^> r (T3 + to) + (t 3 - to) + Z^^Or- - Dy^T)^) • 

r 

(30) 

I 



Here the a matrices are Pauli matrices operating on the spin indices and the r matrices are Pauli matrices operat- 
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ing on the particle-hole indices (tq is the identity matrix 
in the particle-hole space). For the sake of clarity of no- 
tations, we have suppressed the sub-lattice index. 

We now re-write the Fourier transform of the above 
Hamiltonian in a Nambu form to get 

ff MF = c+EE 4A«^> (3i) 

k a, (3 

where C is defined by Eq. 29 and we have now used the 
4-component spinors 



fk,A,/3 fk,B,/3 f- 



■k,A,/3 



f- 



k,B,(3 



(32) 



A, B refer to the two sublattices of the honeycomb lattice 
(as shown in figure 2), j3 =t, | denotes the spin and Hk, a /3 
is given by 



(|-tf)fi a /s Z*p(k) o r a/3 (fc) 

[^a(fc)]* -I> a (-fc) 

o -[Y af} (-k)\* -$-s)n aP -K^(-A)]* 

[I> a (fc)]* -^a(-A) -(|-*)fia/5 

(33) 

We have used the following notations 

£ a p(k) = ~(1 - 2S)e aP (k) - A8e a p(k\ (34) 



For a given spin liquid ansatz, we diagonalize the ma- 
trix Hk as pkDkp\, define the vector 7^ = p\otk and 
determine the values of the mean field parameters 
Dp 1 and rrij. 

To obtain the self consistent solution, we begin with 
an ansatz consistent with magnetic ordering and with 
a combination of hopping and pairing decouplings, and 
allow the system to evolve to a fixed point by self- 
consistent iteration on the values of the mean field 
parameters 28 . As all the mean field parameters are 
quadratic in the fermionic variables, each step in the 
iteration process requires an evaluation of the expecta- 
tion values of quadratic fermion operators in the ground 
state, which are re-calculated iteratively to obtain the 
self-consistent solution. The details of the self-consistent 
equations are given in Appendix A. 

This brings us to the spin liquid ansatz which we de- 
scribe next. 



A. The Spin Liquid Ansatz 

In general we have a nineteen-parameter mean field 
model which needs to be solved self-consistently. These 
fields are: 

On p - links : E* p , D* p , (35) 



(where p = x,y,z) and the on-site magnetization m^. 
In the spin liquid regime, the magnetization is zero and 
we have eighteen complex parameters and the magneti- 
zation. A self consistent mean-field analysis in terms of 
this eighteen(+ one) parameter model suggests that the 
stable mean-field states that we find involve only nine pa- 
rameters or their gauge equivalent forms, in addition to 
magnetization in one phase. Thus we study this nine (+ 
magnetization) parameter model which captures both the 
spin liquid and the magnetically ordered ground states. 
The numerical calculations can further be simplified by 
a correct choice of gauge. To this end, we use insights 
from the exact solution of the Kitaev model. 4 This, as 
shown below, can be obtained by choosing the following 
form for the nine parameters: 

On p — links : Df , Ef G Imaginary, (36) 
DV p G Real, (37) 

(p = x, y, z) with the remaining components set to zero. 
In this gauge, at the Kitaev limit, the dispersion is diag- 
onal in terms of Majorana fermion modes 17 ' 19 as is found 
in the exact solution of the Kitaev model. 4 We use the 
same basis as used by You et al., 17 in which the four 
Majorana fermions are defined as follows 

X = -L(/ it + /^); x l = J- {M -fl) 

X? = ^C/4 + /i); X? = ^(/it"/i)- (38) 

With this ansatz, we now move on to describe the 
two different phases and the phase transition separat- 
ing them. However, before attempting to describe the 
general mean field results, we wish to elaborate on the 
Kitaev limit and the structure of the gauge theory in the 
next two sub- sections. 



1. The Kitaev Limit 

In the Kitaev limit of the model, the above ansatz 
recovers the exact result. 4 Most of the end results in 
this limit are similar to those obtained in Refs. 19 and 
17, because in this limit all these are equal to the ex- 
act solution 4 . However, we point also point out some 
technical differences with our present spinon decomposi- 
tion scheme. The Hamiltonian is given in terms of these 
Majorana fermions as: 



= \ E E ((i - <WK P (x° Xi +P - xlxl +P 

i p 

— (1 — fip,y)iDi jP (XiXi+p ~ XiXi+p + XiXi+p ~ XiXi+p) 

(1 — 3p,x)Di p(Xi Xi+p ~r~ Xi Xi+p ~ X% Xi+p ~ X% Xi+p)^ 

(39) 
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We can rewrite the single occupancy constraint for the 
complex fermions (eq. 10) in terms of the Majorana 
fermions 17 as 

XiXhiXi = \- (40) 

Using this, we can rewrite the spins in terms of the Ma- 
jorana fermions as 

s? = ixhl, sf = i X hl S! = i X iX l (4i) 

which is the original formulation used by Kitaev in the 
solution of his model, with our Majorana fermions nor- 
malized such that {xfiXj} = $ij$ a ^ A set of plaquette 
operators, 

W p = 2 6 5f S% S£Sf SiSl, (42) 

are defined on the individual plaquettes of the lattice, 
where the sites 1 — 6 traverse a honeycomb plaquette as 
shown in figure 2. (The factor of 2 6 which is present 
in our definition of W p is due to the plaquette operator 
being written in terms of spins, rather than Pauli ma- 
trices as in the original formulation of Kitaev. 4 ) These 
plaquette operators commute with the original Kitaev 
spin Hamiltonian and with one another, which allows the 
Hilbert space to be split into eigen-spaces of these oper- 
ators, enabling the exact solution. These operators do 
not commute with the mean field Hamiltonian; however, 
that these operators take the same value in the mean-field 
solution as in the exact solution. 19 

To make a connection with Kitaev's original solution 
we now express our results in terms of the Majorana 
fermions. By construction, the Majorana fermions in- 
troduced in Eq. 38 are the modes in which the band 
structure is diagonal. While x\ X 2 an d X 3 form the flat 
bands, the single dispersing band is made up of the x° 
fermions. 17 In terms of the original solution of Kitaev, the 
dispersing fermion is the single gapless Majorana mode, 
while the flat band fermions describe the frozen Z^ fluxes, 
as we now show. The flat bands arise from the fact that 
the mean-field Hamiltonians for x^X 2 an d X 3 become 
disjoint, i.e., the hopping for these fermions are non-zero 
only on x, y or z bonds respectively. For the hopping on 
the z-link, we have, 

Z(ixh 3 j - iX 3 jXi) (43) 

where 5 is expressed in terms of the mean field param- 
eters and ij are neighbours on a z-link. The eigenvalues 
are given by ±|H|, independent of fc, and therefore these 
form the flat bands. At half filling, the lower energy state 
(lower flat band) is occupied. To compare with the exact 
solution, the Majorana bilinear XiXj nas to be identi- 
fied with the Z2 gauge fields defined on the z-links, u^. 4 
Indeed, we identify 

utj = 2i X h 3 j = *(x?X 3 - X 3 X?)- (44) 



In the ground state, clearly the eigenvalues of u\- are ±1. 
Similarly we can introduce ufj and w(- on x and y links 
respectively. Now we can re-write the flux operators W p 
in Eq. 42 (using 41, 40 and the fact that xfxf = \) as 

W p = 2*SZS$SZSZS%S$ = ul 2 u^ul,ul b ut Q ul x (45) 

It is now clear that in the ground state the plaquette 
operators W p have an expectation value of +1. For a 
small departure from this Kitaev point, one can still use 
the variables vF- and W p . However, these are no longer 
static, but acquire dynamics as the corresponding Majo- 
rana fermions starts dispersing. 

The fermionic mean-field theory of this state describes 
a Z2 spin liquid, as we will show explicitly in the next 
subsection. At the mean field saddle-point, the values of 
different parameters are given by 

- iDV x = E* x = Df >y = Ef <y = Df >x = -iDf z = 0.190608i, 
D* x = -iDV y = E? x = -0.0593918i, (46) 

values which have been determined by self-consistent 
iteration 28 , as described above. 

The resultant spinon spectrum is given in Figure 4. 
There are 8 bands which, characteristic of Bogoliubov 
Hamiltonians, are symmetric about zero energy. The flat 
bands are threefold degenerate. At half filling for the 
spinons the lower four bands (red) are filled while the 
upper four bands (blue) are empty. While the flat bands 
are gapped, the two dispersing bands meet at the bound- 
ary of the hexagonal Brillouin zone with a characteristic 
Dirac spectrum. Hence the spin liquid that we are de- 
scribing is indeed gapless and matches with the spinon 
spectrum obtained in the exact solution of the Kitaev 
model. This provides a useful check on the validity of 
our mean field solution, as well as a controlled limit from 
which we can perturb the model. 

The presence of the pairing term indicates that, in 
terms of the complex fermions, the spin liquid is a "su- 
perconductor" for the spinons. We can analyze the 
symmetry of the pairing amplitude. In order to de- 
termine the properties of the pairing around the Dirac 
node, we isolate the dispersing band by examining the 
X° fermionic modes and returning to the original basis 
of Dirac fermions. For the x° modes, the Hamiltonian is 
given by 

i p 

= f E £(At + 4)( w + /kt) 

i p 

= I E £( M (/^/-fctB + h tA fl w )e-^B + h.c.) 

k p 

(48) 

where M = 0.38122i From here we can expand the 
pairing terms about the K-points in the brillouin zone. 
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FIG. 4: (Color online) Spinon spectrum in the Kitaev limit. 
Bands shown in red are occupied, bands shown in blue are 
unoccupied. 



Defining k f = k + 



2tt \ 



B. The gauge structure 

At this point, before actually discussing the results of 
our mean-field calculations, we wish to discuss the gauge 
structure of the our spin liquid ansatz. 

At the outset, it should be noted that the projec- 
tive symmetry group (PSG) classification including the 
triplet decoupling channels has not been comprehensively 
studied. While a comprehensive discussion of these PSGs 
is beyond the scope of our present work, we indicate the 
relevant issues in the context of the Heisenberg-Kitaev 
model by extending the formalism introduced by Shin- 
dou and Momoi 27 . 

While the formulation outlined above is more suited 
to calculations of the mean field spectrum and self- 
consistent solutions, to decipher the nature of the spin 
liquid and the gauge transformations we wish to cast the 
above decoupling within an SU(2) formalism. 

In order to examine the nature of the spin liquid state, 
it is worthwhile to formulate this Hamiltonian in another 
basis. The transformation into this basis is defined by 



(51) 



M . 



^■dispersing ) g (1 H~ ^ ^ 

MV3 
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( k' x + ik f y ). 



(49) 



However, we would like to emphasize that the above 
chiral p-wave pairing does not necessarily imply time- 
reversal symmetry breaking, which is now implemented 
project ively. 1 The structure of the pairing terms differs 
from the work of Burnell and Nayak 19 , who found p y 
pairing about the Dirac points, by choosing a different 
basis for the fermions which is related to the present one 
by a gauge transformation. 

We can further calculate the spin-spin correlation func- 
tions within mean field theory. Using the Majorana rep- 
resentation we find that this is given by: 



(srsf) 



(Xi Xi ) (Xi X7 / 



(50) 



Since the second correlation function involves absolutely 
flat bands, it is only non-zero when a = f) and when 
i = j or i and j belong to the same unit cell. Hence the 
spin correlation are short ranged even if the spin liquid 
is gapless. This is a novel feature of the Kitaev spin 
liquid, where exact calculations 7 also indicate that such 
correlations vanish beyond nearest neighbour. 

We would like to point out here that, when the model 
is perturbed with the Heisenberg term, the gapped flat 
bands acquire a weak dispersion, but still remain gapped. 
Within perturbation theory, this is expected to lead 
to exponentially decaying spin-spin correlation decaying 
with a length-scale characteristic of the energy-gap. 22 



where the transformation matrix is given by 



A = 



1000- 

1 
10 
0-10 



(52) 



and fi is given by Eq. 28. In the new basis, the f- are 
given by 



fi 



n 



fit fhi fti f 



(53) 



In this basis, we can write the set of gauge transfor- 
mations which leave our physical spin degrees of freedom 
invariant in a block diagonal form, 



Wi 



Vi 
Vi 



(54) 



where the Vi matrices form a two dimensional represen- 
tation of SU(2). The spinon Hamiltonian (Eq 31), when 
written in the new basis, is invariant under the simulta- 
neous gauge transformation 



fi -+ wJi, u' ltP -+ w i+p ul P w}. 



(55) 



where U- p = AUi^A^ gives the analog of Bogoliubov-de- 
Gennes Hamiltonian in the new basis. 

In order to study the low energy degrees of freedom 
in this theory, we allow gauge fluctuations of the U- p 
matrices of the form 



(56) 
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where these k, 1 matrices are block diagonal four by four 
matrices 



k = 







rf 



(57) 



where the rf are Pauli matrices which act on the gauge 
degree of freedom, and generate our gauge transforma- 
tions. We also take note of a set of matrices which gen- 
erate our spin rotational symmetry, which in this basis 
are given by 



'I, 



(58) 



where the a 1 are again Pauli matrices, acting on the spin 
degrees of freedom of our fermions. 

To determine the gauge structure, we now consider the 
product of the U' ip matrices around different lattice loops 
based at any site i, 1 



(59) 



c 



Here the product is taken over the bonds of the loop 
Ci, beginning and ending at the site i. Using the given 



notations, we can write our U' ip matrices as 



U- £ 



(60) 



where E a = (E°, E 1 , E 2 , E 3 ), ^ = («°, k\ k 2 , k z ) (E° 
and k° are 4x4 identity matrices and the other ma- 
trices are given by Eqs. 57 and 58), and the are 
complex numbers. We note that, unlike the singlet case, 
in the triplet decoupling scheme both the gauge and the 
spin generators enter in tJ[ p . (In other words, in the 
singlet decoupling 1 E° is the only spin generator that 
enters since the decoupling channels are invariant under 
spin rotations). 

Similarly the loop function can be written as 



(61) 



where the are determined by the values of the 
on the links of the loop Ci. 

If all the loop functions, based at any site, can be 
rotated into a form which commutes with one or more 
gauge generators, then the set of such generators form the 
invariant gauge group (IGG) and the low energy gauge 
fluctuations belong to the IGG. 1 

Taking our ansatz in the Kitaev limit of the model, the 
structure of the U' ip matrices is given by 

tJ[ = -i^ 3 E 3 + iQk 2 Y? + iVE 1 , (62) 



The loop functions in this limit (in our choice of gauge) 
have a typical structure which is given by 

P(d) = T(^°E° - /^E 1 + k 2 Y? + ^ 3 E 3 ), (64) 

where T is a constant. We find that these cannot be 
brought into a form which commutes with any of the 
gauge generators, and hence the only kind of low energy 
gauge fluctuation allowed has the form Wi = e l6iK where 
€i = or 7r. This gives an IGG = Z2 and we have a Z2 spin 
liquid. This spin liquid has gapless Major ana excitations 
(see previous sub-section) and is indeed the Kitaev spin 
liquid. This IGG survives throughout the entire regime 
of S in which magnetic order is absent. This completes 
our discussion regarding the invariant gauge group of the 
gapless spin liquid, which we have now shown to be Z2, 
as expected from Kitaev's exact solution. 4 

However, in the presence of magnetic ordering the 
above picture is no longer true. Magnetic order drives 
the E? p parameter to zero, which changes these U' i p ma- 
trices into a form given by 

^ = zQ^ 2 E 2 +P^ 1 E 1 , (65) 

where 



p = k- 

8 lv 2 



5) + 4(5(1 

8) + 45{l-6 p , v )]DV f 

The loop functions are now given by 

P(d) = T(k°S° + k 3 S 3 ), 



(66) 
(67) 

(68) 



where f is a different constant. It is now easy to show 
that the ft 3 matrix commutes with these products, and 
the gauge transformations of the form Wi = e lBiK (6i G 
[0, 27r)) leave our ansatz invariant. Thus, in this case, 
the IGG= U(l) and hence the low energy gauge fluc- 
tuations are described by a compact U{1) gauge theory 
which has a gapless photon and also instantons (space- 
time monopoles) where the gauge flux may change in 
integral multiples of 2tt. In addition, using the gauge 
transformation 



fL 



ia y 
ia y 



fl 



i,B-> 











Kp ( 69 ) 



the Dij can now be completely rotated into the vec- 
tors (when Ei iP is zero) hence explicitly showing that this 
is a U(l) spin liquid. Later, we shall see that the spinon 
spectrum is gapped (in the presence of magnetic order- 
ing), and that this state is actually a gapped U(l) spin 
liquid which is unstable to confinement in two spatial 
dimensions. This significance of this instability will be 
discussed later. 



where 



Q=W-5p,y)Dl r 



(63) 



IV. THE RESULTS OF THE MEAN FIELD 
THEORY AND BEYOND 

We now discuss the results of our mean field calcula- 
tions as a function of 5. These mean field results are 
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FIG. 5: (Color online) (a) The spinon band structure for the Heisenberg-Kitaev model with nonzero Heisenberg coupling and 
zero net magnetization, at the point 5 — 0.3. (b) The spinon band structure for the Heisenberg-Kitaev model with nonzero 
Heisenberg coupling and non-zero net magnetization, at the point 6 = 0.23. 



obtained by solving the self-consistent equations (given 
in Appendix- A) for the various mean field parameters. 

a. Region surrounding the Kitaev limit (5 « 0.5): 
Near the Kitaev limit, we find that the spinon bands 
which were flat in the pure Kitaev limit gain a dispersion, 
with energy which scales with the distance from the ex- 
actly solvable point. These bands do not contribute sig- 
nificantly to the low energy theory due to the fact that 
they remain fully gapped, and the low energy spinon ex- 
citations remain consistent with those of the pure Kitaev 
model. Figure 5(a) shows the band structure in this re- 
gion, at the value of S = 0.3. As the strength of the 
Heisenberg coupling is increased, the mean field parame- 
ters show only a slight change prior to the onset of mag- 
netic ordering. 

b. Region with non-zero magnetic order (5 < 0.26): 
At S ~ 0.26, we find that the system begins to admit 

a non-zero magnetic order parameter as a self-consistent 
solution. The magnetic order parameter jumps discon- 
tinuously to a finite value at this point, indicating a first 
order phase transition. The spinon band structure differs 
significantly in this phase, which includes the formation 
of a band gap. Figure 5(b) shows a cut of the spinon 
bands in this region, at the value of S = 0.23. We also 
note that the values of all of the mean field parameters 
are changed by this ordering, and that the value of the 
hopping order parameters are driven to zero. As we con- 
tinue to increase the strength of the Heisenberg coupling 
(decrease S) we see that the magnetic order parameter is 
increased, and the pairing amplitudes are driven to zero 
as well (below 5 « 0.15). Once the pairing amplitudes 
are zero, all the spinon bands become flat (not shown) 
and also have an energy gap. 

The self-consistent values of the different mean field 
parameters are plotted in figure 6. This shows that the 
magnetic order parameter discontinuously turns on at 



0.8 r— 
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FIG. 6: (Color online) The magnitude of the mean field pa- 
rameters, plotted as a function of 6. The onset of magnetic 
order triggers a first order phase transition. The symbols are 
a guide for the eye. 



0.26. Below this value of S there is finite magnetic 



order. Also at that value of r5, 
tinuously. 



Ef- goes to zero discon- 



A. Interpretation of mean field results 



As discussed, for S > 0.26 (i.e. a > 0.76) there is 
no magnetic order and the E{ iP and Di jP fields are non- 
zero. This, as we have already discussed already, is a 
Z2 spin liquid which is continuously connected to the 
exactly solvable kitaev spin liquid (obtained for a = 1). 
This has gapless Majorana fermion excitations and short 
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range spin correlations. 

On the other hand, for 5 < 0.26 there is magnetic or- 
der. However, at the mean field level spinon bands are 
well defined and there are dispersing spinons. Only the 
D i p fields are non-zero in this phase. However, we have 
already shown (see Eq. 69) that these D iyP fields can be 
gauge rotated into E ijP fields and hence this regime repre- 
sents a U(l) spin liquid. Since the spinon band structure 
is gapped, we essentially have a gapped U(l) spin liquid 
with ferromagnetic order. We can call this phase a FM* 
in order to distinguish it with the regular ferromagnetic 
order. In addition to the gapped fermionic spinon exci- 
tations, there is a gapless emergent U(l) gauge photon 
present in this phase. This arises from the underlying 
U(l) gauge fluctuations that this phase allows. However, 
contrary to our regular electrodynamics, this emergent 
U(l) gauge group is compact 1 . Hence, as noted earlier, 
in addition to the photon, it allows instanton processes 
where the magnetic flux changes by integral multiple of 
2tt. This turns out to be significant, as we discuss in 
the next subsection. We also point out that in general 
the FM* phase has a Goldstone mode (spin wave) that 
is in general gapped because the spin-rotation symme- 
try is broken explicitly (except at the points 6 = and 
(5 = 0.5). 

Since within mean field theory magnetic order turns on 
discontinuously, the transition is first order within the 
limits of our numerical resolution. The jump is about 
20% of the saturation value. Hence, within mean field 
theory, we have a first order transition between the Z^ 
spin liquid with gapless spinons with Dirac dispersion to 
a FM* phase with gapped spinons. 



B. Beyond mean field theory: Instantons and 
confinement of FM* 

In this sub-section, we discuss the issue of instability 
of the FM* phase to a conventional ferromagnetic phase. 
As already pointed out, the compact U(l) gauge the- 
ory that describes the low energy excitations of the FM* 
phase allows tunnelling processes where the magnetic flux 
changes (instanton events). It is known from the work 
of Polyakov 20 that in (2 + 1) dimensional compact U(l) 
gauge theory, when the matter fields carrying electric 
charges(spinons) are gapped, the instanton events are al- 
ways relevant. Thus, once we incorporate fluctuations to 
our mean field solutions, we have to take into account the 
effect of such tunnelling processes. Once such instanton 
events are taken into account, the spinons, which carry 
gauge charges, are confined to gauge neutral objects- the 
spins. This confinement is not, however, a straightfor- 
ward consequence of magnetic order, as a stable gapless 
U(l) phase with deconfined spinons and coexisting mag- 
netic ordering can occur in two spatial dimensions 29,30,37 . 
We emphasize that it is the U{1) gauge structure of the 
magnetically ordered phase, combined with the gapped 
nature of the spinon excitations, 20 which is responsible 



for the confinement through the proliferation of instanton 
events. 

The above discussion indicates that once we move be- 
yond mean field theory and take the instantons of the 
compact U(l) gauge theory into account, the spinons in 
the FM* phase undergo confinement. However, the ferro- 
magnetic order parameter would survive due to the fact 
that it is gauge invariant. Such a confined phase is con- 
tinuously connected to the regular ferromagnetic phase 
for the spins and we end up with a FM phase (or the 
stripy phase for the unrotated spins). Thus, we indeed 
get a direct transition from the Z2 spin liquid to a stripy 
phase, albeit discontinuously. 



V. DISCUSSION 

We now summarize our results. We have obtained a 
slave-particle description of the HK model and used it 
to describe the phase transition between a spin liquid 
and the magnetically ordered stripy phase within slave 
particle mean field theory. In the Kitaev limit of the 
model, we have shown that this formulation reproduces 
the expected excitation spectrum and that the plaque- 
tte operators which enable the exact solution are in a 
vortex free configuration in the ground state. Upon the 
inclusion of a small non-zero Heisenberg term we have 
found a similar low energy theory, although the bands 
which are dispersion-free in the Kitaev limit gain a dis- 
persion. We have analyzed the gauge structure of the 
model, and have seen that in the absence of magnetic 
order the Z2 IGG which describes the Kitaev spin liq- 
uid state remains the IGG of our ansatz. The magnet- 
ically ordered phase that we get by destroying the Z 2 
spin liquid is, within mean field theory, a gapped U(l) 
spin liquid which has stripy magnetic order. However, 
existing results imply that such a spin liquid is unstable 
to confinement, which immediately drives a transition 
to the regular stripy antiferromagnetic phase. Within 
mean-field theory, the above transition turns out to be 
discontinuous. Our description allows for a coherent de- 
scription of the spin liquid and the magnetically ordered 
phase as well as of the phase transition connecting them. 

The present numerical results 8,15 cannot conclusively 
shed light on the nature of the above transition between 
the spin liquid and the stripy antiferromagnet. However, 
these results seem to suggest that the transition is ei- 
ther continuous or weakly first order. While our mean 
field theory indicates a first order transition, we are re- 
quired to incorporate quantum fluctuations beyond the 
mean field to address the issue of a possible continu- 
ous transition between the phases. In fact it is some- 
what easy to see why the transition appears to be first 
order in our present calculations. Once we neglect the 
gauge fluctuations, we can treat the fields E^- and 
as "order parameters" along with the actual magnetiza- 
tion order parameter m^. A Landau-Ginzburg theory in 
terms of these fields can then be obtained by integrat- 
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ing out the fermions. Such a "multi-order parameter" 
Landau-Ginzburg theory with repulsive interactions be- 
tween the order-parameter densities generically gives a 
first order transition within mean-field theory. 31 Hence, 
the results of our mean field calculations can be under- 
stood within this framework. However, a shortcoming of 
such a naive Landau-Ginzburg analysis is the fact that it 
cannot take into account the effect of gauge fluctuations. 
Once such fluctuations are accounted for, the E^- and 
Y>ij fields can no longer be treated as order parameters, 
since they are not gauge invariant, and so the above naive 
Landau-Ginzburg theory breaks down. This opens up the 
possibility of subtle gauge fluctuation effects driving this 
transition to second order. It is known from earlier stud- 
ies in frustrated magnets that such "Landau forbidden" 
generic continuum quantum phase transitions may occur 
(e.g. deconfined quantum critical points 32 ) where naive 
mean field considerations break down or do not apply, 
since there are no local order parameters (e.g. Topo- 
logical phase transitions 33 ). Hence such a second order 
transition would be unconventional and potentially in- 
teresting, particularly in context of the the possibility of 
realizing the HK model in material systems. 8 ' 9 Hence, 
this would be an ideal subject of a future study. Similar 
transitions in spin rotation invariant systems have been 
recently studied both numerically 34,35 and from the field 
theoretical perspective. 36 
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Appendix A: Self-consistent Mean Field Theory 

In order to determine the values of our mean field pa- 
rameters, we consider the matrix H^, defined as 



Hh = 



(Al) 



where Hk, a p is given by 33. We diagonalize as 
U/cD/cU^ and define the vector 7^ — U^a/c, where is 
given by 



*fc,t k,i 



(A2) 



where a\ a is given by 32. For a given value of the mean 
field parameters we evaluate these for each wavevector 
and determine the values of the mean field parameters 
E^(p) and D^(p) in this state. To do this, we must look 
at expectation values of the form (fj a fip) and (fj a fip). 
We have assumed previously that the parameters E and 



D are uniform throughout the lattice, and we use this 
fact to rewrite these expectation values as 



7— links 



N, it 



N. 



k 

EEE^™^* ,m,m' 1 k : m' I 



ik-R 



k V m' 



Y Y U k,l,l' U k,m,l> 



k I' (unoccupied) 



(A3) 



for 1 corresponding to /3, A, m corresponding to a, £?, and 
the sum over unoccupied V referring to the eigenvectors 
of Hk corresponding to positive eigenvalues. Similarly, 



-(fi,pfi+-y,a) 



^ (fi,pfi+'y,a) 



j— links 



N, it 



N. 



site 



EEE^w^* 

,m,m Ik, m/ 



t \-ik-R 



k V m' 



E E U k ^U* w e 



k I' (unoccupied) 



(A4) 



for 1 corresponding to /3, A and m corresponding to a, B 
on the -k portion of the vector. We also need to evaluate 
terms of the form (fj a fia) to determine our magnetiza- 
tion m. As above, we determine these by evaluating 



^site 

% 

= y2(fka,Afka,A) 

■lS ftit.p. 



site 



Y( U h,l^il' U k,l,rnHk,rn>) 



k I' m' 



N, 



E E itfc.w- 



k I' (occupied) 



(A5) 



where 1 corresponds to a, A. We can determine the occu- 
pancy of the sites on sublattice B similarily. Using these 
expressions, we can determine the expectation value of 
our mean field parameters beginning with any ansatz. 
Doing this iteratively (with the previous solutions as our 
new ansatz at each step) allows us to determine an ap- 
proximate self-consistent solution to our mean-field equa- 
tions. 
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